function subregionTfcPeaks(subjectinfo, useSelectedElectrodesTf)

if ~exist('subjectinfo', 'var')
    Subjects;
    subjectinfo = subject.CW;
end

if ~exist('useSelectedElectrodesTf', 'var')
    useSelectedElectrodesTf = true;
end

%% Prepare data
if useSelectedElectrodesTf
    Electrodes = subjectinfo.Electrodes(logical([subjectinfo.Electrodes.SelectedManually])); % Use only 'Selected' electrodes
else
    Electrodes = subjectinfo.Electrodes;
end
Electrodes = Electrodes(~[Electrodes.RejectInPreprocessing]);
% uniqueSubregions = unique({Electrodes.Subregion});
uniqueSubregions = unique({Electrodes.Region});
uniqueSubregions = uniqueSubregions(~strcmp(uniqueSubregions, '')); % Throw out the empty string...
for i = 1:length(uniqueSubregions)
    selections(i).subregion = uniqueSubregions{i};
    selections(i).subregion_fieldname = uniqueSubregions{i};
    selections(i).subregion_fieldname(uniqueSubregions{i}=='-') = '_';
%     selections(i).electrodes = ElectrodeInfo.FindBySubregion(Electrodes, uniqueSubregions{i});
    selections(i).electrodes = ElectrodeInfo.FindByRegion(Electrodes, uniqueSubregions{i});
end

% setToDo = {'SphereUnambiguousResponse', 'SphereAmbiguousResponse', 'FaceHouseUnambiguousResponse', 'FaceHouseAmbiguousFaceResponses', 'FaceHouseAmbiguousHouseResponses'}'
setToDo = fieldnames(subjectinfo.EpochData);
for i = 1:length(setToDo)
    epochset = setToDo{i};
    DATA = load(subjectinfo.EpochData.(epochset).TfcDatafile);
    DATA.NrOfEpochs = subjectinfo.EpochData.(epochset).nrOfEpochsFound;
    tfc = BuildTfcSet(DATA, selections);
    disp('Building dataset...');
    [datatable header] = TfcToDatatable(tfc);
    disp('Writing dataset to disk...');
    save([epochset '.mat'], 'datatable', 'header', '-v7.3');
%     xlswrite([epochset '.xls'], [header;datatable]);
end

bins = {'4-7Hz', '7-13Hz', '13-30Hz', '30-70Hz', '70-130Hz'};
pptCfg = plotPeakTiming_defaultCfg();
for b = 1:length(bins)
    figure;
    pptCfg.binSelection = bins{b};
    for i = 1:length(setToDo)
        epochset = setToDo{i};
        load([epochset '.mat'], 'datatable', 'header');
        subplot(length(setToDo),1,i);
        plotPeakTiming(datatable, header, pptCfg);
        title(epochset);
    end
    a = axes('position',[.02 .02 .96 .96],'Visible', 'off');
    text('Position', [0.01 0.99], 'String', [subjectinfo.Name ' / ' bins{b}], ...
        'HorizontalAlignment', 'left', 'VerticalAlignment', 'cap', 'Parent', a);
    text('Position', [0.5 0.01], 'String', datestr(now, 'yyyy-mm-dd HH:MM:SS'), ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'baseline', 'Parent', a);
end;
end

function tfc = BuildTfcSet(D, selections)

tfc.freqbindescriptions = D.tfc.freqbindescriptions;
tfc.freqbins = D.tfc.freqbins;
tfc.times = D.tfc.times;

for i = 1:length(selections) % For every subregion
    disp(selections(i).subregion);
%     disp(length(selections(i).electrodes));
    % Gather data
    for frq = 1:length(D.tfc.freqbindescriptions) % For every frequencybin
        tfc.data.(selections(i).subregion_fieldname)(frq).frequencyrange = D.tfc.freqbins(frq,:);
        if length(selections(i).electrodes)>1
            m = []; s = [];
            for k = 1:length(selections(i).electrodes) % For all selected electrodes
                m(:,k) = D.tfc.data.(selections(i).electrodes(k).Label)(frq).mean;
                s(:,k) = D.tfc.data.(selections(i).electrodes(k).Label)(frq).std;
            end
            tfc.data.(selections(i).subregion_fieldname)(frq).mean = mean(m,2);
            tfc.data.(selections(i).subregion_fieldname)(frq).std = sqrt(sum(s.^2, 2));
        else
            tfc.data.(selections(i).subregion_fieldname)(frq).mean = D.tfc.data.(selections(i).electrodes(1).Label)(frq).mean;
            tfc.data.(selections(i).subregion_fieldname)(frq).std = D.tfc.data.(selections(i).electrodes(1).Label)(frq).std;
        end
        PEAKS = scorePeaks(tfc.data.(selections(i).subregion_fieldname)(frq), tfc.times - 5000, num2str(str2double(D.NrOfEpochs) * length(selections(i).electrodes)));
        tfc.data.(selections(i).subregion_fieldname)(frq).peaklocs = PEAKS.peaklocs;
        tfc.data.(selections(i).subregion_fieldname)(frq).peaktype = PEAKS.peaktype;
        tfc.data.(selections(i).subregion_fieldname)(frq).firstEffect = PEAKS.firstEffect;
        tfc.data.(selections(i).subregion_fieldname)(frq).pvalue = PEAKS.pvalue;
    end
end

end